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Disordered packings of soft grains are fragile mechanical systems that loose rigidity upon lowering 
the external pressure towards zero. At zero pressure, we find that any infinitesimal strain-impulse 
propagates initially as a non- linear solitary wave progressively attenuated by disorder. We demon- 
strate that the particle fluctuations generated by the solitary- wave decay, can be viewed as a granular 
analogue of temperature. Their presence is manifested by two emergent macroscopic properties ab- 
sent in the unperturbed granular packing: a finite pressure that scales with the injected energy 
(akin to a granular temperature) and an anomalous viscosity that arises even when the microscopic 
mechanisms of energy dissipation are negligible. Consistent with the interpretation of this state 
as a fluid-like thermalized state, the shear modulus remains zero. Further, we follow in detail the 
attenuation of the initial solitary wave identifying two distinct regimes : an initial exponential decay, 
followed by a longer power law decay and suggest simple models to explain these two regimes. 



I. INTRODUCTION 

The defining property of solid materials is their abil- 
ity to resist external mechanical perturbations, as em- 
bodied in their finite elastic moduli. However, it is the 
geometry and topology of a materials' architecture that 
ultimately controls the strength of the elastic response. 
For example, if we remodel a solid block of steel into 
a granular aggregate of steel balls just in contact with 
their nearest neighbours, both the linear clastic moduli 
and sound speed of the aggregate drop to zero. These 
extreme softness, that we term fragility, stems directly 
from the anharmonic interaction between macroscopic 
grains in contact. It is independent of the material the 
grain is made of which could be as hard as steel or as 
soft as bubbles. Weakly connected harmonic networks of 
cross-linked polymers or covalent glasses can also become 
fragile, irrespective of the individual bond strength, pro- 
vided that the number of mechanical constraints is too 
low to maintain rigidity. 

The critical point at which a packing of soft repulsive 
grains become fragile (unjams) can be accessed experi- 
mentally by progressively reducing the pressure, or al- 
ternatively the average overlap, to zero. Right at the 
critical point, even the tiniest mechanical strains must 
propagate in the form of supersonic solitons and shocks 
since the linear sound speed is zero. Nesterenko coined 
the term sonic vacuum to designate this state in a gran- 
ular chain and demonstrated the existence of solitonic 
excitations both experimentally and theoretically [l|, y] . 
The fate of these strongly non-linear excitations in amor- 
phous two or three dimensional granular media near the 
(un)jamming transition or in weakly connected polymer 
networks is only gradually emerging Q. 

The very existence of these non-linear excitations and 
the associated critical points impinges on a very precise 



tuning of some geometrical parameters like the particle 
overlap that must vanish at random close-packing or the 
coordination number of the network that must be tuned 
exactly to maintain rigidity. It is natural to enquire 
how fluctuations (that we expect to be violent given the 
absence of linear restoring forces) around these special 
points change their critical behavior. However, the study 
of fragile nano-mechanical systems affected by quantum 
fluctuations is still in its infancy. Even the impact of 
thermal fluctuations on the un-jamming transition is not 
very well studied, partly because granular media, that 
represent the key experimental arena for these investiga- 
tions, are clearly athermal. 

In this article, we study numerically and analytically 
the decay of solitary wave excitations in two dimensional 
amorphous or mass-disordered packings of grains that are 
just in contact with their nearest neighbors, see Fig. ([T]a- 
b). For weak mass disorder, the solitary wave excitation 
generated in response to an impulse, decays exponen- 
tially at early times, with a rate that depends upon the 
amount of disorder. In the long time limit, the initially 
well defined solitary wave soon transitions into a trian- 
gular shock like profile, decaying as a power law, with 
a universal exponent consistent with ^ and independent 
of the amount of disorder. This power law decay is the 
dominant mechanism of attenuation in hexagonal lattices 
with strong mass disorder as well as in jammed packings. 
Finally, we study the fluctuating state (similar to a ther- 
malized state) that emerges after the complete disinte- 
gration of the initial solitary wave, see the inset of Fig. 
([T|d) and the trail of the soliton in Fig. ([TJ:). The resulting 
fluid-like state has two emergent macroscopic properties 
absent in the unperturbed granular packing: (i) a flnite 
pressure that scales with the energy originally injected in 
the form of solitary wave and (ii) a finite viscosity, that 
arises despite no microscopic dissipative mechanism are 
present in our simulations. 
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FIG. 1. (a) A schematic illustration of how a hexagonal pack- 
ing at zero pressure (with vanishing longitudinal and trans- 
verse sound speeds) and weak mass disorder, responds to an 
impact at one of its ends, by generating a solitary wave exci- 
tation, (b) As the solitary wave propagates and interacts with 
the weak mass inhomogeneity, its amplitude decays. The de- 
cay at early times is nearly exponential and the process of en- 
ergy loss is such, that it excites several smaller solitary waves. 
Over the long run, a noisy state (velocity fluctuations) span- 
ning the size of the packing emerges (zoomed region), (c) The 
initial solitary wave is now no longer identifiable. Instead we 
observe a triangular shock like profile that is defined as an 
envelope over the smaller excitations. In this regime, the de- 
cay of the leading edge follows a power law. Eventually, no 
leading propagating edge is visible. 



II. ATTENUATION OF NON-LINEAR 
SOLITARY WAVES 



To study the dynamical response near the critical 
point, we perform molecular dynamics simulations for 
N = 4000 soft frictionless particles arranged in a hexag- 
onal lattice or as a jammed amorphous packing. The par- 
ticles interact pair-wise with the Hertz potential V{6) ~ 

^ (tr) ' '^h^^'2 K sets the energy scale and the par- 
ticle overlap S is measured in units of the average disk 
radii R. By discretizing the packings into bins along the 
X— direction (defined as the longitudinal direction here), 
we impart an initial speed Up to the left most bin and 
obtain the subsequent dynamics by integrating the equa- 
tions of motion using the velocity Verlet method |15l |. 



We follow how the initial energy imparted to the pack- 
ing evolves into a wave profile that propagates along the 
direction of initial impact (see schematic Fig. ((!])) by av- 
eraging out the transverse speeds over the bins. 

In Fig. ([2]), we show the time evolution of the solitary 
wave excitation generated in response to an impulse for 
three cases: (a) a hexagonal packing with grain masses 
distributed as a normal random variable with a small 
variance e (b) the same hexagonal packing with large 
variance in the masses of the grains and (c) a jammed 
amorphous packing prepared at a vanishingly small pres- 
sure P~ 10~^. In all three cases, we find that in response 
to an impulse, the mechanical excitation initially evolves 
into a well defined solitary wave with a spatial extent 
around 5 grain diameters, that propagates along the di- 
rection of impact. Remarkably, the solitary wave excita- 
tion is to a good approximation, the soliton-likc solution 
first seen by Nesterenko in a one-dimensional granular 
chain of macroscopic non-cohesive particles interacting 
via the Hertz-potential 043] ■ 



A. Exponential decay 

As the solitary wave propagates through a weakly dis- 
ordered hexagonal packing (see inset of Fig. ([2^), it be- 
gins to attenuate and the initial stages of this attenua- 
tion is well approximated as an exponential decay. The 
isotropic elasticity of the hexagonal packing allows us 
to model the solitary wave as a one dimensional exci- 
tation - we make use of the quasi-particle approxima- 
tion to understand its decay in this regime |14{ . Con- 
sequently, we model the interaction of the solitary wave 
quasi-particle with the mass inhomogeneity as an elas- 
tic collision process (conserving momentum and energy) 
that results in the disintegration of the solitary wave into 
two new solitary waves during each collision. After prop- 
agating through n grains (undergoing multiple collisions) 
we find that the mean energy of the solitary wave is given 
by (see App. A) 



Tn 



(1) 



where, Tq is the initial energy of the solitary wave, Tn 
is the solitary wave energy after traversing n beads and 
e is the variance in the masses of the beads, assumed to 
be distributed as a Gaussian random variable. As shown 
in the inset to Fig. ([Da), we find this estimate to be 
in very good agreement with our numerical observations 
on an hexagonal packing and weakly-disordered granular 
chains M. 



B. Power law attenuation. 

By contrast, for the hexagonal packing with strong 
mass disorder Fig. ([2]b), the exponential regime is barely 
identifiable and the solitary wave begins to evolve into a 



shock like triangular profile. For the decay in this regime, 
we find a striking similarity in all three cases : the long 
time decay in a hexagonal packing with weak mass disor- 
der, the dominant regime of decay in a hexagonal packing 
with strong mass disorder (large e) as well as the dom- 
inant mechanism of decay in amorphous jammed pack- 
ings, see Fig. ([5]a-c). In all cases, a triangular shock like 
profile emerges, whose leading edge decays as a power law 
a;"'" with an exponent approximately r « 0.5 (solid red 
line) . This exponent is thus independent of the amount 
of disorder. 

Notice, the shock like profile is obtained as an envelope 
over several smaller excitations that span the system up 
to the leading propagating edge. Thus, unlike the initial 
solitary wave excitation whose spatial extent is around 
5 grain diameters, the shock solution spans several hun- 
dreds grains. In a recent study, the existence of such 
long wavelength triangular shock like solutions in a one 
dimensional chain of granular particles was established 
[201 ■ Here, we find that due to material inhomogeneities, 
an initial excitation in two dimensional disordered pack- 
ings, naturally evolve into long wavelength shock solu- 
tions. As we outline below, the power law decay follows 
as a simple consequence of energy conservation. 

In the long wavelength approximation, the average 
strain field 5 satisfies the continuity equation 5t+5^''^5x = 
0, where subscripts denote partial derivatives with re- 
spect to time t and space x [2^. This equation has a 

similarity solution 5{x,t) ^ (j) . Upon integrating once 
with respect to x and differentiating once with respect 
to i, we obtain the corresponding solution for the parti- 
cle velocity field (j){x,t) ~ (j) . Since the energy of the 
beads enclosed within the shock envelope is conserved, it 

follows that E ^ L^ dx<j>'^{x,t) ~ jrw = contant, where 
Xf is the position of the shock front and consequently. 
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yielding the exponent i. As shown in Fig. (I2]a-c), solid 
red lines, we do find good agreement with the numeri- 
cally determined exponent r « 0.5 for the decay of the 
shock front in the power law regime. 



III. EMERGENT HYDRODYNAMIC STATE 

The simple model for the exponential decay of the soli- 
tary wave suggests that at each subsequent collision with 
an inhomogenity, the solitary wave splits approximately 
into two solitary waves- a leading pulse (the main degree 
of freedom that is damped exponentially as a result) and 
a smaller solitary wave, that is either reflected or trans- 
mitted depending upon the mass ratio. Therefore, once 
we allow sufficient time for the leading solitary wave to 
disintegrate completely such that a leading pulse is no 
longer distinguishable from the background, we expect to 
reach a state comprised of several smaller solitary waves 



with different energies. The solitary waves in turn inter- 
act with each other in-elastically, (as a consequence of 
the Nesterenko equation of motion being non-integrable 
[21) and thus their interaction may be thought of as inher- 
ently dissipative. However, in addition to these processes 
that seek to distribute the energy initially concentrated 
in a single solitary wave excitation into multiple smaller 
solitary waves, the structural disorder in two dimensional 
amorphous packings also spills a part of the energy into 
transverse degrees of motion. Through a series of such 
intrinsic dissipative mechanisms, we eventually reach a 
fluctuating equilibrium-like state that spans the entire 
flnitc-sizcd packings under-investigation. 

In order to rationalize the physics behind this fluctu- 
ating state that at first appears to be just noise, we as- 
sume that there is a well defined long wavelength, small 
frequency collective hydrodynamical mode with wave- 
number k in the x-direction. Upon binning in the x 
direction, we define the coarse grained particle current 
density j(r,t) = ^W X]i=i ■^i(0'^('"^i"j(^)) and its Fourier 
transform j„(k,t) = -^J2tl^^o.ity^"'^'^■ (The bold 
face is a vector notation used for spatial coordinates 
whose cartesian components are denoted by a). Upon 
setting k ~ {k, 0) and a = x or y, we determine nu- 
merically the corresponding longitudinal and transverse 
current density auto-correlation functions as Ci^tik,t) = 
(j;*t(k, 0)j;,t(k, i)), where the angular brackets denote en- 
semble averaging over the initial time. The longitudi- 
nal and transverse power spectral densities Pi^t{k,co) are 
evaluated numerically using fast Fourier transform from 
the respective current density auto-correlation functions 
as 



Pi.tik,^) 



dt 



'Ci,t{k,t). 



(2) 



In Fig. ([3| b) , we plot the longitudinal (red squares) 
and transverse (red circles) dispersion curves obtained 
numerically from the power spectral density Fig. ([3[ a) 
in the emergent fluctuating state. For comparison, we 
plot the corresponding longitudinal (black squares) and 
transverse (black circles) dispersion curves for a highly 
compressed jammed packing far from the critical den- 
sity, prepared at a pressure of P~ 10~^. Since the total 
potential energy of a jammed packing is related to its 
pressure via E~ P^'"^, we choose an initial impact speed 
Up w 2.0 to generate a solitary wave in the weakly com- 
pressed packings (P~ 10^®). This leads to a total initial 
solitary wave energy Esw that is comparable to the en- 
ergy of the highly compressed packing Epc and thus al- 
lows for a more meaningful comparison between the two 
states. 

As seen in Fig. ([3|b), highly compressed jammed pack- 
ings behave as ordinary solids with a finite bulk and shear 
modulus and this translates into a finite sound speed 
manifest in the linear regime of the longitudinal (black 
squares) and shear (black circles) dispersion curves. In 
contrast, exciting a jammed packing prepared at a very 



small pressures (that a-priori has zero longitudinal and 
shear sound speeds), leads to a linear dispersion regime 
for the longitudinal modes, but no such regime is ob- 
tained for the transverse modes. The slope of the linear 
regime in the longitudinal dispersion curves, corresponds 
to the speed of long wavelength hydrodynamical sound 
modes. Defining the sound speed c as the second deriva- 
tive of the induced potential energy; Esw'i leads to the 
relation, c ~ ^sw UM' closely matching the numerical 
data in Inset to Fig. (I3]b), red squares. Thus, the speed of 
the long wavelength hydrodynamical sound modes scales 
with the energy of the initial solitary wave Esw injected 
into the system. 

This is analogous to the scaling relation for pre- 
compressed jammed packings at a finite packing frac- 
tion S, where the sound speed scales as c ~ eJq , with 
Epc being the potential energy due to the finite pre- 
compression 6. Thus, in so far as longitudinal sound 
modes are concerned, a rigidity induced by statically 
compressing a marginally compressed packing is analo- 
gous to the rigidity induced by exciting a marginally com- 
pressed packing with a finite energy wave. Note there- 
fore, one can easily replace the source of energy by a 
heat bath and thereby obtain a thermally induced rigid- 
ity upon making the substitution E— >■ ksT. However, 
unlike a state that is truly in thermal equilibrium, an 
external perturbation over the fluctuating state created 
by the disintegration of a solitary wave, will further raise 
its energy due to the absence of a fluctuation-dissipation 
mechanism. The emergent state is thus at best described 
as a quasi-equilibrium state. 



In contrast to the longitudinal modes, the transverse 
modes obtained by energizing a marginally compressed 
packing do not show a well defined linear regime, see 
Fig. ([3]b), red circles. This is in stark contrast from 
a statically compressed jammed packing where a linear 
transverse dispersion regime (owing to the finite shear 
modulus) with a slope that scales with the amount of 
pre-compression (and does not depend upon the solitary 

wave energy injected) as c '^ Ep(-, is expected (Fig. ([3] 
b), black circles). Thus, the shear modes excited by in- 
jecting energy into a packing near its critical point are 
purely diffusive and the medium does not develop a finite 
shear modulus. There is therefore a profound difference 
between the resultant states obtained by (a) statically 
compressing a granular packing near its critical point ver- 
sus (b) injecting energy either in the form of an excitation 
or thermally. In the former case, one obtains a solid-like 
medium with finite bulk and shear moduli, while in the 
latter, we ffuidize the medium. This means that, even af- 
ter the non-linear wave excitations (characteristic of the 
jamming point) are strongly attenuated, one cannot de- 
scribe the system purely in terms of the normal modes 
of a (linear elastic) solid because the system has de facto 
become a fluid. This is one of the key conclusions of our 



work and it heralds the signal property of packings at the 
jamming threshold: they are fragile. 

In Fig. ([3] c) , we show how the half width (inverse life- 
time r ) of the longitudinal hydrodynamical modes de- 
pend upon the wavenumber k. We find that it scales 
anomalously as t~^ ^^ fc^'^. For purely hydrodynami- 
cal modes obeying the Navier-Stokes equation, the half 
width scales with the wave number as hw ~ 77 A;^, where 
77 is the shear viscosity. However, extensive numerical 
and analytical studies have shown that in one dimension 
the time correlation functions (whose long time integrals 
by definition correspond to macroscopic transport prop- 
erties such as diffusivity and shear viscosity) do not de- 
cay exponentially but display long time tails, decaying 
as power laws instead |16l4l8| . This phenomenon indi- 
cates the breakdown in low dimensions of the standard 
continuum assumption embedded in the Navier-Stokes 
equation - the strength of fluctuations is too strong for 
simple coarse-grained theories to hold. 

Analytic studies based on mode coupling theory pre- 
dict that for d — 1 dimensions, the shear viscosity itself 
scales with the wave number (in the small wave num- 
ber limit) as r]{k) ~ fc"^/'^ and thus to leading order, we 
obtain a half width that scales as hw-^ rj{k)k'^ = k^/^. 
These results have also been validated numerically in 
more recent studies on one dimensional non-linear spring 
chains [IllliJI- The predicted 5/3 power law dependence 
is consistent with our findings in Fig. (|3]c) solid red line. 
Our packings are two-dimensional but the emergent hy- 
drodynamic description is effectively one-dimensional be- 
cause we are injecting compressional solitary wave fronts 
in the x directions only, see Fig. (1). Thus, despite 
the absence of any microscopic dissipativc mechanism in 
our simulations, there is an emergent shear viscosity that 
scales with the wave number and gives the longitudinal 
and shear hydrodynamical modes a finite lifetime. 



IV. CONCLUSION 

To conclude, we find that in response to an impulse, 
hexagonal packings with mass disorder and amorphous 
packings, initially generate a well defined solitary wave 
excitation that is progressively attenuated. For weakly 
disordered hexagonal packings, the amplitude of the soli- 
tary wave excitation decays exponentially at early times, 
with a rate that depends upon the amount of disorder. In 
the long time limit, we observe a transition to a power law 
regime, with a universal exponent that is independent of 
the amount of disorder and is the dominant mechanism 
of decay in the strongly disordered case, where an expo- 
nential decay is no longer identifiable. We find that the 
physics of the solitary wave attenuation in jammed amor- 
phous packings at their critical point corresponds to the 
strongly disordered regime in hexagonal packings where 
an initially well defined solitary wave soon transitions 
into a triangular shock like profile, decaying as a power 
law. We understand the two stages of decay using simple 



analytical models and scaling arguments. We then study 
the resultant fluctuating state (similar to a thermalized 
state) that emerges after the complete disintegration of 
the initial solitary wave. This fluid-like state has emer- 
gent macroscopic properties absent in the unperturbed 
granular packing: (i) a finite pressure that scales with 
the energy originally injected in the form of a solitary 
wave and (ii) a finite viscosity, despite no microscopic 
mechanism of dissipation being present. 
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Appendix A: The exponential decay 

Consider a lattice of beads of mass m, where nearest 
neighbours interact with a one sided non-linear poten- 
tial of the form V{S) = —6°', with S being the average 
overlap between neighbouring beads, k the force constant 
and a > 2 the exponent of the non-linear potential. If 
the initial overlap; 5 — >■ 0, the effective spring constant 
defined as the second derivative of the potential vanishes, 
and as result, the lattice does not support linear sound at 
zero temperature. Consequently, mechanical strains gen- 
erated in response to an impulse evolve into long-lived 
non-linear solitary waves and the solitary wave energy E 

and momentum P satisfy the classical relation E = t^ , 

where moff ~ l.3m is the mass of the solitary wave quasi- 



particle as a function of the mass of the beads to [l|, Q ■ 
Here, the momentum P = ttIcsU where, U physically 
corresponds to the solitary wave amplitude. 

The quasi-particle approximation of the solitary wave 
has been used to study the disintegration of a solitary 
wave across a mass interface [a,[ij|. Consider an interface 
between two regions of sonic vacuum with grain masses 
77ii,TO2 respectively. For mass ratios A — ^^^ close to 
1, a solitary wave initially moving with amplitude Uq is 
seen to split into two new solitary waves, whose ampli- 
tudes f/i./, f/2./ niay be obtained using an elastic collision 
model that conserves the quasi-particle energy and mo- 
mentum - 



rrii^cffUQ = mi^offUij + TO2,cffC^2,/j 



from which one can easily derive 



In order to study the propagation of the solitary 
wave in a medium with weak mass inhomogeneity, 
consider a chain of beads with the mass ratio of neigh- 
bouring beads i,j related via rrii = Aijmj, where 
Aij = 1 + Nij(0,e^), The normal random variable 
Nij(0,e^) has mean and variance e^. Upon appealing 
to the localized nature of the solitary wave (its width 
being around 5 bead diameters and also independent 
of the amplitude of the solitary wave), we treat each 
bead as an interface and invoke the quasi-particle elastic 
collision model. 

Thus, the energy of the leading solitary wave after the 
first collision is 



Ti = 



l + eNi(0,l) 
(l + fNi(0,l))- 



:Tn 



1 + eNi(0, !))(! - eNi(0, 1) + -e^Nf (0, 1) 



Tn 



l--e'NliO,l) 



To, 



where, we have retained terms up to order e^. Iterating 
this process n— times, we find that the solitary wave 
energy after it has propagated n beads diameters is 



T„ 



l-ieX(0,l))--(l-J^'N?(0,l) 



Tn 



Retaining terms only to order e^, we find 



To 



fc=i 



l--/x'in), 



where, X^i^) is the chi function with an expectation value 
n. Taking the mean, we obtain that the average solitary 
wave energy after propagating n beads diameters reads 



To 



A - ine2 



(A3) 



This result agrees quantitatively with numerical simula- 
tions in ID [Sl and qualitatively with our observations in 
2D packings. 
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FIG. 2. (a): The time evolution in response to an impulse in 
a hexagonal packing with weak mass disorder e. An impulse 
response generates a solitary wave excitation that decays ex- 
ponentially at early times. In the long time limit, a shock 
like triangular profile emerges whose leading edge decays as 
a power law with an exponent x'"'^ (solid red line). (b)The 
time evolution in response to an impulse in a hexagonal pack- 
ing with strong mass disorder e. An impulse response gener- 
ates a solitary wave excitation but the regime of exponential 
attenuation is barely identifiable. Instead, a shock like trian- 
gular profile soon emerges whose leading edge again decays 
as a power law with an exponent a;~°'^ (solid red line), (c) 
The time evolution in response to an impulse in a jammed 
amorphous packing prepared at a pressure of P~ 10~® and 
an average overlap between grains So ~ 10~^. An impulse re- 
sponse still generates a solitary wave excitation but like (b), a 
shock like triangular profile soon emerges whose leading edge 
decays as a power law with an exponent a;~°* (solid red line). 




FIG. 3. (a) The power spectral density for longitudinal modes 
that emerge in marginally compressed amorphous packings 
{P ~ 10~®) by the complete disintegration of an initial soli- 
tary wave excitation, (b) The red squares shows the numer- 
ical data for the longitudinal dispersion curves. The slope 
of the linear regime scales with the energy of the solitary 
wave as c ~ _E^' ^^ for Hertzian interaction, that we identify 
with long wavelength longitudinal hydrodynamical modes. In 
contrast, the red circles shows the numerically obtained data 
for the transverse dispersion curve where no linear regime 
is seen. The shear mode is therefore non-propagating. For 
comparison, shown are linear dispersion curves for longitu- 
dinal (black squares) and transverse (black circles) obtained 
for highly compressed jammed packings, prepared at a pres- 
sure P~ 10~^.(c) The half-width obtained numerically from 
the longitudinal modes as a function of wavenumber on a 
linear scale and compared against the analytical estimates 
h.w.~ fc^" (solid red line). 



